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ABSTRACT 

We present multidimensional simulations of the early convective phase preced- 
ing ignition in a Type I X-ray burst using the low Mach number hydrodynamics 
code, MAESTRO. A low Mach number approach is necessary in order to perform 
long-time integration required to study such phenomena. Using MAESTRO, we are 
able to capture the expansion of the atmosphere due to large-scale heating while 
capturing local compressibility effects such as those due to reactions and ther- 
mal diffusion. We also discuss the preparation of one-dimensional initial models 
and the subsequent mapping into our multidimensional framework. Our method 
of initial model generation differs from that used in previous multidimensional 
studies, which evolved a system through multiple bursts in one dimension before 
mapping onto a multidimensional grid. In our multidimensional simulations, we 
find that the resolution necessary to properly resolve the burning layer is an order 
of magnitude greater than that used in the earlier studies mentioned above. We 
characterize the convective patterns that form and discuss their resulting influ- 
ence on the state of the convective region, which is important in modeling the 
outburst itself. 

Subject headings: convection — hydrodynamics — methods: numerical — stars: neutron- 
X-rays: bursts 

1. Introduction 

Type I X-ray bursts are possibly the most frequent thermonuclear explosions in the 
universe and provide a large amount of observational data that can be used to determine 
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the properties of matter near the surface of a neutron star. To make meaningful inferences 
about these properties from observational data, however, we mu st have a proper theoretical 
understanding of the bursting phenomena (jBhattacharyyall2010l ). The basic XRB paradigm 
takes place in a mass-transferring, low-mass X-ray binary (LMXB) system in which the 
neutron star's companion has filled its Roche lobe and is dumping H- and/or He-rich material 
onto the surface of the neutron star. Depending on the acc retion rat e and composition, there 
are several burning regimes that will trigger an XRB (see iBildstenl (120001 ) for an overview) . 
The general idea is that a column of accreted material — or heavier-element ash from prior 
stable burning of accreted material — builds up until the temperature sensitivity of the energy 
generation rate at the base of the layer exceeds that of the local cooling rate and a thin-shell 
thermal instability forms. The instability eventually causes a runaway of unstable burning 
resulting in an outburst. 

One-dimensional hydrodynamic studies reproduce many of the observable features of 
XRBs such as burst energies (~ 10 39 erg), rise times (seconds), durations (IP's - 100's of 



seconds) and recur rence ti mes (hours to days) (IWooslev fc Weaver! 11984 iTaam et al.lll993 
Heger et al. 2007a ; also see Strohmayer fc Bildsten ( 20061 ) for a review of XRBs). By con 



struction, however, one-dimensional models assume that the fuel is burned uniformly over 
the surface of the star, which is highly unlikely given the la rge disparity between the ther- 
malization and burning timescales of the accreted material (ISharalll982l ). Furthermore, the 
Rossi X-ray Timing Explorer satellite has observ ed coherent oscillations in the lightcurves 
of ~ 20 outbursts from LMXB systems (first by IStrohmayer et al.l Il996l ; more recently by 
Altamirano et al.ll2010l and references therein). The asymptotic evolution of the frequency of 



such oscillations suggests they are modulated by the neutron star spin frequency (IMuno et al. 



20021 ). Oscillations observed during the rising portion of an outburst lightcurve are therefore 
indicative of a spreading burning front being brought in and out of view by stellar rotation. 
Additionally, oscillations observed during the decay phase of the burst are thought to be 
caused by unstable surface modes that may depend critically on the local heating and cool- 
ing rates during the burst (INarayan fc Cooperll2007l . and references therein). The manner in 
which the burning front spreads and propagates throughout the accreted atmosphere is not 
well known, and a proper multi dimensional modeling of the conditions in the atmosphere 



prior to outburst is needed (e.g. iFryxell fc Woosleylll982bl ). 



Prior to the actual outburst, the burning at the base of the ignition column drives con- 
vection throughout the overlying layers and determines the state of the material in which 
the burning front will propagate. One- dimensional simulations of XRBs usually attempt to 
parameterize the con vective overturn and mixing using astrop hysical mixing-len gth theory 
(jBohm- Vitensd 1 1 9581 ) or through various diffusive processes (see iHeger et al.ll2000l for a thor- 



ough discussion). Recent multidimensional simulations of stellar convection (see lArnett et al. 



-3 - 



20091 . and references therein), however, show a large discrepancy in, for example, the velocity 
of a typical convective eddy when compared to one-dimensional models in the case of stellar 
evolution codes that use mixing-length theory. Indeed, there has recently been an effort put 
forth in the astrophysical community, the so-called Convection Algorithms Based on Simu- 
lations, or CABS, to derive from multidimensional simu lations a more phys ically motivated 
prescription for handling convection in one dimension ( lArnett et al.ll2008l ). To date, such 
methods have not propagated into the XRB-simulation community and a proper treatment 
of convection, without assumptions, requires simulation in multiple dimensions. 

Multidimensional simulations of any aspect of XRBs, however, have hitherto been rather 
restrictive. A burning front can propagate either supersonically as a detonation or sub- 
sonically as a deflagration . Full hydrodynamic XRB detonation models in the spirit of 
Fryxell fc Woosleyl ( Il982al ) or lZingale et al.l ( 1200 ll ) require a thick (~ 100 m) accreted helium 
layer. Such deep layers are only produced by very low accretion rates, which are inconsistent 
with the majority of rates inferred from observations of XRBs, and therefore the burning 
front in most XRBs likely propagates as a deflagration. Deflagration models are difficult to 
compute with standard compressible hydrodynamics codes due to the long integration times 
required. One possible solution is to eliminate the effect of acoustic waves in the system, 
allowing the time step to be controlled by the fluid velocity, rather than the sound speed. 
Such a method can be derived using low Mach number asympt otics; classic examples o f low 
Mach numb er approaches incl ude the incompressible, anelastic f Qgura fc Phillips Il962[) and 
Boussinesq flBoussinesql Il903[ ) approximations. To this end, ISpitkovsky et al.l ( 120021 ) used 
a simple, shallow-water, 2-layer, incompressible fluid to model the vertical structure of a 
deflagration front and showed how rotation coupled with convection may play an important 
role in regulating the spread of the front over the surface of the neutron star. 



More recently, iLin et al.l (120061 ) developed and applied a low Mach number approxima- 
tion method to the problem of convective burning at the base of an accreted layer in an 
XRB system. Their method, however, was first order accurate in space and time and did not 
allow for the evolution of the hydrostatic base state, a feature that is needed to capture the 
expansion of the atmosphere in response to heating. Furthermore, Lin et al. did not model 
the surface of the accreted layer, which is vital to understanding bursts that exhibit pho- 
tospheric radius expansion (PRE bursts); such bursts are crucial in determining the stellar 
properties of neutron stars ( Steiner et al. 201ol . and references therein). 

In this study we use MAESTRO ( Nonaka et al.|[201ol ). a multidimensional low Mach number 
hydrodynamics algorithm for astrophysical flows, to model the convection leading up to an 
outburst of a pure 4 He accretor. MAESTRO is a second-order accurate, conservative method 
that uses rectangular grid cells. The algorithm is capable of capturing the expansion of the 
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atmosphere due to large-scale heating, while capturing local compressibility effects due to 
reactions, thermal diffusion and compositional changes. A semi-analytic method is used to 
generate one- dimensional initial models. These models are then augmented by being evolved 
in a one-dimensional stellar evolution code; this evolution allows for the approximation 
of the convective cooling and leads to a model that is closer to satisfying the thin shell 
instability condition. The resulting model is then mapped into MAESTRO and evolved in 
multiple dimensions. 

The main goals of this paper are to explore and describe the challenges of modeling 
XRBs in multiple dimensions and to better understand the convective phase that precedes 
the outburst. The remainder of this paper is as follows. In Section [2J we describe the gener- 
ation of the initial models and the subsequent mapping into a multidimensional framework. 
In Section [31 we describe the MAESTRO algorithm, including the addition of two new modules 
not present in the original algorithm. Specifically, we have added thermal conduction and a 
"volume discrepancy" correction term to the velocity field to ensure that the solution does 
not diverge from the equation of state. In Section HJ we describe the results of our multi- 
dimensional simulations. In particular we discuss the resolution needed to properly resolve 
the burning layer, the effects of including the thermal diffusion and volume discrepancy cor- 
rection terms, the expansion of the base state due to heating and finally the nature of the 
convective behavior and its effect on the atmosphere. We conclude in Section [5] by summa- 
rizing our results and describing our plans to extend the current study to mixed H/He XRB 
sources. 



Initial Models 



We begin our calculations by generating a one-dimensional initial model of the accreted 
layer in hydrostatic (HSE) and thermal equilibria on the surface of a neutron star. We assume 
a plane-parallel geometry — that is, the gravitational acceleration, g, is assumed constant 
throughout the domain, which is justified because the thickness of the accreted layer (~ 10 
m) is much less than the radius of the neutron star (~ 10 km). We assume a 4 He layer is 
accreted on top of a 56 Fe neutron star with a trace abundance (10 -10 ) of 12 C. We choose 
a pure 4 He accretor both because the corresponding nuclear reaction network, 3a burning, 
is simple compared to the slow, /3-decay-limited burning processes in bursts involving H, 
and beca use ultra-compa ct XRB sources are possible pure 4 He accretors (4U 1820-30, for 



example; ICumming 



2003m. We include the forward and reverse 3a re action rates as g i yen in 



Caughlan fc Fowler! (119881 ) wit h electron screening contribu tions from lGraboske et al.l (119731 ) 
for the weak regime and from lAlastuey fc Jancovicil (119781 ) for the strong regime. 
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There are several approaches to one-dimensional model generation in the literature. In 
our approach, we begin with a semi-analytic initial model and then augment the model 
to account for convective cooling. We also discuss proper mapping of the one- dimensional 
model to our multidimensional framework. 



2.1. Semi-Analytic Models 

The semi-analytic approach to model generation involves integration of the heat equa- 
tion and an entropy equation, 

at 1 q,^t? 

(1) 
(2) 

where c is the speed of light, a the radiation constant, k the opacity (including radiative 
and conductive contributions), T the temperature, F the o utward heat flux and dy = —pdr 



dT 


3kF 


dy 


AacT 3 


dF 






= o, 


dy 



with y(r) the column-depth (see lCumming fc Bildsten!l2000l for details of this method). Note 



that ([T]) can give a thermal profile that is superadiabatic — in practice, the thermal gradient 
is restricted to be dT/dy < (dT/dy) s where the subscript s means along an adiabat. Also 
note that for simplicity, equation fl2]) neglects any compressional heating contributions from 
the accretion itself and assumes the accreted material is not burning during the accretion 
phase — this is a steady-state configuration. There is, however, an outward heat flux from 
pycnonuclear reactions deep within the neutron star crust; we approximate this flux as a con- 
stant value throughout the accreted layer, F = 200 keV per nucleon. The integration starts 
at the top of the 4 He atmosphere (arbitrarily at y top = 10 3 g cm -2 ) where a radiative zero 



solut ion is assumed, and continues until the thin shell instability condition (IFushiki fc Lamb 



19871) 



de^a de coo \ 

It > ~dW' [6) 

is reached at y = ?/base- The local cooling rate is typically approximated from ([T|) and fl2]) as 

^cool ~ T, n- l^J 

6ny z 

When ([3]) is attained, the composition for y > ybasc is switched to 56 Fe and integration of 
([I]) and (j2D resumes until a thick enough substrate is formed such that y base is sufficiently 
far from the bottom of the computational domain, y{r = 0) = 10 12 g cm -2 in our studies. 

The approximation, (j4]), works well in one dimension because the only efficient way 
the system can cool (neglecting weak reactions) is via conduction and radiation, which 
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enter through the opacity. When more spatial dimensions are added to the system and 
there is heating from below from nuclear reactions, the fluid is free to overturn and cool 
via convection. Now we have a situation where the local multidimensional cooling rate, 
e cooi, muiti-d = e cooi+ e conv, exceeds the initial approximation and ([3]) may no longer be satisfied. 
Therefore, such a semi-analytic model is no longer close to runaway and to evolve the system 
in multiple dimensions until ([3]) is reached is intractable even with the advantages of a low 
Mach number approximation code. 



2.2. Kepler-supplemented Models 

One way to overcome the difficulties with evolving the model described in the previous 
section in multiple dimensions is to explicitly include an effective convective cooling term 
in the approximation to the local cooling given by equation (j4j). This effective convective 
cooling can be included via mixing-length theory typically found in stellar evolution codes. 
Using the semi-analytic model described above a s initial conditions, the one- dimensional 



stellar evolution code, Kepler (I Weaver et al .111 9781 ). was used to construct the remainder of 



the underlying neutron star with R ns = 10 km and M ns = 1.87M© (Woosley, 2010; private 
communication). The system is then allowed to evolve in one dimension whereupon nuclear 
burning heats the base of the layer, and the convection prescription develops a well-mixed 
and nearly adiabatic region of 12 C ash overlying the 4 He base. This results in a model that 
is much closer to satisfying the thermal instability criterion, fl3]), when mapped into multiple 
dimensions. 



2.3. Mapping to Multiple Dimensions 



The data from Kepler are given in a Lagrangian (mass) coordinate system and we need 
to convert them to an Eulerian (ph ysical) coordinate sy stem for use in MAESTRO. We use 
a procedure similar to that found in IZingale et al.l (120021 ) to ensure our initial model is in 
HSE. Given the density, temperature and composition from the Kepler evolution, we call 
the equation of state to get the pressure. We then discretize the HSE equation and solve for 
the non-uniform Eulerian grid spacing corresponding to the Lagrangian grid points, 



Ti-l 



1 pi- pi-i 



9 V2 (Pi + Pi-l) 



(5) 



where r is the radial coordinate, p the pressure and p the density. We set r± — to complete 
the description of the grid. The transition from the pure 56 Fe neutron star (at r tra ns) to the 
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4 He atmosphere (at r trans+ i) is a step function as a result of the initial Lagrangian data. Such 
sharp transitions can be a source of numerical noise and oscillations as the solution evolves 
on an Eulerian grid. To minimize the numerical noise, we smooth the interface by adding n 
uniformly distributed coordinate points between r trans and r trans+ i. The temperature at these 
new points is linearly interpolated between T trans and T trans+ i. Then X( 4 He) and X( 12 C) at 
the new points are filled with a tanh profile: 



4>i = a tanh + C (6) 

V <P J 

where a = (0trans+i - Rtrans) /2, r c = (r trans + r trans+ i) /2, </> c = (0 trans + 0t™ns+i) /2 and y? is 
a parameter to set the smoothness. X( 56 Fe) is then found from the constraint ^2 k Xk = 1, 
and p and p are found by using an iterative Newton-Raphson technique with the equation of 
state and §5§ at these new points. This smoothed model is then linearly interpolated onto a 
completely uniform grid, with rj = r^-i + Ar, and is again put into HSE using fl5]) and the 
equation of state. Values of n = 50 and <p = 3 were used to smooth the models presented in 
this work. 

Figure [T] shows the result of this procedure for two models that were evolved in Kepler 
until the base of the 4 He atmosphere had reached a temperature of 3.67 x 10 8 K (solid line, 
hereafter referred to as the cold model) and 5.39 x 10 8 K (dotted line, hereafter referred to 
as the hot model). The density at the base of the 4 He layer for the cold model is 1.4 x 10 6 g 



cm 3 and is 1.2 x 10 6 g cm 3 for the hot model. For comparison, the initial model of lLin et al. 
( 120061 ) had a base temperature and density of 2 x 10 8 K and 4 x 10 6 g cm -3 , respectively 



The cold model has a peak in 12 C production around r = 382 cm (i.e., the base of the 4 He 
layer in both models) that appears smoothed in the more evolved hot model. Both models, 
however, have an extended region of well-mixed 12 C that extends to r = 624 cm (r = 812 
cm) for the cold (hot) model. These initial models contain no multidimensional velocity 
information from the Kepler simulations. We therefore make no assumptions about the 
nature of the convection when the models are mapped into multiple dimensions in MAESTRO. 



3. Hydrodynamics Algorithm 

For our multidimensional simulations, we use the low Mach number stellar hydrody- 
namics algorithm, MAESTRO. This code is appropriate for flows in the low Mach number 
regime, where the characteristic fluid velocity is small compared to the speed of sound. Note 



valid for such flows. A series of papers (see 



Almgren et al.l (l2006bl ) — henceforth Paper II, 



Almeren et al. 




2006a) 


Almeren et al. 


J2008 





-henceforth Paper III, 
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and IZingale et al.l ( 120091 ) — henceforth Paper IV) describe the derivation of the low Mach 
number equation set, its algorithmic implementation, and the initial application to con- 
vection in a white dwarf preceding a Type la supernova. We use the most recent version 



of the algorithm, which includes local adaptive mesh refinement, as given in iNonaka et al. 



(|2010|)— henceforth Paper V. 



One key advantage of using a low Mach number approach is the increase of allowable time 
step size, which enables long-time integration. Standard compressib le hydrody namics codes 



for a strophysical applications such as CASTRO (lAlmgren et al. 



20101 ) or FLASH flFrvxell et al. 



20001 ) evolve a fully compressible equation set, i.e., the Euler equations, which allows for 
the formation and propagation of shocks. For low speed convective motion in our pre-burst 
convection studies, we do not need to explicitly follow the propagation of sound waves. Our 
low Mach number equation set does not contain acoustic waves, and therefore MAESTRO is 
able to take time steps constrained by the maximum fluid velocity, rather than the maximum 
sound speed. As an example, if the maximum Mach number of the flow is M ~ 0.01, 
we will obtain a factor of 1/M ~ 100 increase in time step size compared to a standard 
compressible approach. Another advantage of a low Mach number method is that the overall 
HSE of the state can be guaranteed by the inclusion of a base state in HSE in the low Mach 
number equation set. This removes the difficulties of maintaining HSE commonly found in 
compressible hydrodynamics codes. 

MAESTRO solves a system of advection-reaction-diffusion equations with the equation 
of state formulated as an elliptic constraint on the velocity. MAESTRO uses a higher-order 
Godunov method to discretize the advective terms, Strang-splitting to couple the reaction 
terms to the advective terms, and a semi-implicit treatment of the diffusion terms. The 
diffusion term and the divergence constraint are formulated as linear systems which are 
solved iteratively using multigrid. The evolution of the one- dimensional base state density is 
also computed. The base state density represents the average state of the atmosphere, and 
is coupled to the base state pressure via HSE. The base state density has its own evolution 
equation that computes the expansion of the atmosphere due to heating and is discretized 
using a higher-order Godunov method. We note that MAESTRO is second-order accurate in 
space and time. 



3.1. MAESTRO Details 

We now provide additional details of the low Mach number equation set and numerical 
implementation. The interested reader is referred to Papers I-V for full details. We use 
a two-dimensional Cartesian formulation with x the horizontal coordinate and r the radial 
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coordinate. The low Mach number equation set is 
d( P X k ) 



dt 



-V ■ (pX k U) + pCj k , (7) 



^ = -U-VU-^-t^ger, (8) 
at p p 

^ = -V-(phU) + ^ + pH mc + V-(k th VT), (9) 

where U, h and fc th are the velocity, specific enthalpy, and thermal conductivity, respec- 
tively. The species are represented by their mass fractions, X k , along with their associated 
production rates, uj k , and H nuc is the nuclear energy generation rate per unit mass. Using 
low Mach number asymptotics (see Paper I) the total pressure, p(x,r,t), is decomposed into 
a base state pressure, p Q (r,t), and a perturbational, or dynamic, pressure, ir(x,r,t), such 
that 1 7r | / p = 0(M 2 ). The base state density, po(r, t), is in HSE with the base state pressure 
such that Vpo = ~Po9 e r, where e r is the unit vector in the outward radial direction. 

Thermal conduction was not present in Paper V, so we have developed a semi-implicit 
discretization for this term. We include full algorithmic implementation details in Appendix 
IA1 and a verification test problem in Appendix IA.1I 

Mathematically, this system must still be closed by the equation of state, which is 
expressed as a divergence constraint on the velocity field (see Paper III), 

V • = ft f S - =MsrY (10) 



r lPo dt 



where (3 is a density-like variable, 



Po(r,t) =p(0,t)exp Qf dr'^j , (11) 

and Ti(r) is the average of Ti = (dlnp/dlnp) s where the subscript s means the derivative 
is taken at constant entropy. We will use an overline notation to represent the average of 
a quantity, which computationally is the arithmetic average of all grid cells at a particular 
radius. The expansion term, S, in ( TTUj) accounts for local compressibility effects resulting 
from nuclear burning, compositional changes, and thermal conduction: 

S = aH nuc + -aY,^k + —Y,P^ + -V-(k th VT), (12) 
h PPp h P 



where £ k = (dh/dX k ) pT(x .^ k) , p p = (dp/dp) TXk , p Xk = (dp/dX k ) Tp(X] j ^ k) and a = 
Pr/ipCpPp) with p T = (dp/dT) pXk , and c p = (dh/dT) pXk . 
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Another addition to the MAESTRO algorithm is the use of a "volume discrepancy" cor- 
rection. Because (fT0|) is a linearization of the nonlinear constraint imposed by the equation 
of state, the thermodynamic pre ssure, Peos — p(p, h, X k ), may drift from the base state 



pressure, po, fjPember et al.lll998f ). To correct for this drift, ( ITU]) is augmented with a term 



that drives the thermodynamic pressure back to the that of the base state: 

'■<* D i-*( , -^-isr^)' (13) 

where / is the volume discrepancy correction factor and < / < 1. In Section l4~4l we explore 
the effectiveness of this term at keeping the overall solution in thermodynamic equilibrium. 

To track the evolution of the base state density, we first define the expansion velocity 
as the average outward velocity: 



w (r,t) = (U-e r ). (14) 

As described in Paper V, we compute Wq by integrating a one-dimensional divergence con- 
straint, found by taking the average of equation ( {TBI : 

d {PqWq) = ^f-g ^_^Po f Po- Peos \ ^ 



Or 



The evolution equation for the base state density can be found by considering the average 
of the continuity equation: 

^ = -V • (p w e r ), (16) 

which we discretize with a higher-order Godunov method. 

We use special care in dealing with the low density region of our simulation. The density 
spans many orders of magnitude, and due to conservation of momentum we may generate 
large velocities in the upper atmosphere that do not affect the solution in the higher-density 
region. Unfortunately, these large velocities reduce the efficiency of our method by reducing 
the time step size. The first technique we use to address this problem is the use of a cutoff 
density, p cu toff, which is the value we hold the density to outside the star. The second 
technique we use is the use of an anelastic cutoff density, p a neiastic, below which we determine 
/So by keeping the ratio Po/Po constant in the divergence constraint in order to minimize 
spurious wave generation. Full implementation details for the cutoff densities are described 
in Appendix A. 5 of Paper V. In this paper, we use p cu toff = Paneiastic = 10 4 g/cm 3 . 

The third technique adopted for the low density region is sponging (or damping), which 
is used to reduce gravity waves at the surface of the star. This technique is commonly 
used in the atmospheric modeling community as lateral boundary conditions of limited area 
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simulations (see, for example, iKesel fc Winninghofflll972t iPerkey fc Kreitzberpl976f ) as well 
as upper boundary conditions to re duce wave reflection off of sharp gradients in the atmo - 
spheric structure (see, for example, iDurran Sz Klemplll983l ; |Purranlll990l ; IChen et al.ll2005l ). 
In addition, we have previously utilized the sponging technique in the study of convection in 
the cores of white dwarfs (Paper IV). Full details for the sponge implementation in MAESTRO 
can be found in Papers III and IV, but in summary, we add a forcing term to the velocity, 
which effectively damps the velocity so that U ncw — > JJ ncw * /damp- We use the following 
formulation for the sponge: 



1. 



r < r 



damp 



damp, min 



COS 



7T 



r tp - 



1 + /, 



damp, min J ; 



/damp,] 



spi 

sp < r < r tp, 

r tp < r, 



r. 



;i7) 



where in our simulations r sp is the radius at which po = 25p cuto ff, r tp is the radius at which 

PO = PcutofT, and /damp, min 

0.010 In Figure [TJ the vertical grey lines correspond to the 
location of r sp and the vertical black lines correspond to the location of r tp for each of the 
initial models. Figure [2] shows the initial profile of the sponge for the cold model. Note that 
as the system evolves it is free to expand thus changing the location of the density cutoffs 
and consequently the location and extent of the sponge. The inclusion of a sponge layer 
does not strictly conserve kinetic energy in the sponged region. We note, however, that the 
material above the surface of the star is at relatively low density compared to the material 
in the convective region, and therefore the total amount of energy non-conservation is small. 
Furthermore, as shown in Figure 4 of Paper III, the inclusion of a sponge layer in the low 
density region of a simulation does not affect the dynamics of the flow in the convective 
region of interest. 



4. Results 

We describe below the results of mapping the Kepler-supplemented models into MAESTRO 
in two dimensions, (x, r), and the system's subsequent evolution. Section fl~Tl describes the 
resolution requirements needed to properly resolve the burning layer. In Section U2] we show 
how the inclusion of thermal diffusion affects the nuclear burning layer and its location. 
We show in Section 14.31 how utilizing a time-dependent base state allows us to capture the 
expansion of the atmosphere due to heating. Section 14.41 shows how including a volume 



1 Note that the form of this sponge is similar to that presented in Section 4.3.1 of Paper III but with 
nAt = 1 at each time step. 
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discrepancy correction keeps the base state thermally consistent with the equation of state. 
Finally we discuss the extent and evolution of the convective region in Section 14.51 

To map the one-dimensional model into MAESTRO, we copy laterally across the domain 
such that (j)(x,r,t = 0) = O ne-d( r ) for each variable <fi in the model. In the following 
analysis, the subscript "max" refers to the maximum value of a quantity in the computational 
domain at a given time step. In two dimensions, we define the average as a function of 
radius, (0) = (</>)(r, t), of a quantity (f> by 



1 N 



ill, 
'J 



where 0™- = <fr(xi, rj, t m ) and N is the total number of grid zones in the lateral, x, direction 
at height Tj at time t m . 



We use the general equation of state of iTimmes &: Swestyl (120001 ) . which includes con 



tributions from electrons, ions, and radiation. We calculate opacities using Frank Timmes' 
publicly available conductivity routine, wh ich includes con tributions from radiation and elec- 



tron conduction processes as explained in iTimmesi (120001 ) . It is important to note that the 



method for calculating opacities used in MAESTRO is not the same as what is used the Kepler 
code; it is likely, however, that the different metho ds give opacities that agree to within a 



factor of ~ 2 (see discussion in iHeger et al.ll2007bl ). The boundary conditions for all simu- 



lations are periodic in the ^-direction to mimic a laterally extended convection region. The 
upper r boundary is outflow to allow for free expansion of the atmosphere. The lower r 
boundary of dense neutron star material is set to a wall with no normal flow. To solve the 
thermal diffusion contribution at the upper and lower boundaries we impose the Neumann 
condition dh/dn = 0, where n is the outward facing normal vector; the enthalpy boundary 
conditions are periodic in the lateral directions. We note that the upper and lower domain 
boundaries are sufficiently far from the burning layer so that they do not affect the dynamics 
of the convection. An advective CFL number of 0.7 was used in all of our simulations. 

As previously mentioned, we do not obtain any multidimensional velocity information 
from the Kepler models; our system is initially static. For convection to begin, the symmetry 
of the system must, therefore, be broken. This can be accomplished either by placing a small 
perturbation at the base of the 4 He layer or by allowing numerical noise from the multigrid 
solver to seed the convective cells. For the simulations presented here, neither approach 
is advantageous over the other, both giving quantitatively similar steady-state convective 
flow fields; we utilize both approaches in our studies and when perturbing we place a small 
(AT/T = 10~ 5 ) Gaussian temperature perturbation laterally centered at height r = 384 cm 
to break the initial symmetry of the problem. 
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4.1. Resolution Requirements 



To date, the only othe r paper in the l iterature regarding multidimensional simulations 
of XRBs as deflagrations fjLin et al.l 120061 ) used a finest resolution of 5 cm zone -1 . They 
presented multidimensional results at 5, 7.5, and 10 cm zone -1 resolutions and remarked 
that there is a "tendency toward convergence with increasing resolution" based on the time 
to reach the peak energy generation rate. It is important to note that our initial models 
are different from those of iLin et al.l (120061 ). In particular, their models only considered two 
species — the accreted layer was pure 4 He and the underlying neutron star was composed 
entirely of 12 C. This caused their models to have a smaller jump in mean molecular weight 
across the neutron star/accreted layer boundary compared to our models. Furthermore, the 
initial conditions for their multidimensional studies were from the results of a one- dimensional 
diffusional-thermal code that evolved the system thr ough sever a l burs ts. These differences 
from our method of initial model generation give the ILin et al.l (120061 ) models an extended 
(~ 100 cm) thermal peak compared to our narrow (~ 10 cm) peak (compare our Figured] 
to their Figure 2). 

The burning layer at the base of the accreted material is very thin; high resolution is 
required to properly model this region. The peak of the thermal profile for the hot model 
is broader than the corresponding peak in the cold model. Consequently, the burning layer 
in the hot model is thicker than that of the cold model — we therefore focus our study of 
resolution requirements on the more restrictive of the two, the cold model. The top panel 
of Figure [3] shows the (H nuc ) profile at t = 1 m s for simulations of the cold model using 



the same resolutions as in the ILin et al.l ( 120061 ) study. Even at this early time there is a 
25% spread in the peak value of (H nuc ) for these resolutions. The bottom panel shows the 
same profile but at several higher resolutions. The peak value of (H nnc ) for the 4 cm zone -1 
simulation is comparable to the peak values in the top panel, but we only see numerical 
convergence of the peak value as we go to higher spatial resolution. In addition, the shape 
of the profile near peak converges with increasing resolution; the 0.25 and 0.5 cm zone -1 
resolution simulations look qualitatively similar. We therefore claim that the burning layer 
is not properly resolved in our models unless a resolution of 0.5 cm zone -1 is used. It is 
important to note that even though our initial models differ, this resolution requirement is 
an order of magnitude higher than what has been previously presented in the literature and 
therefore significantly increases the computational cost of our XRB simulations. 

Under-resolving the burning layer artificially boosts the energy generation rate, which 
in turn over-drives convection. Figure H] shows a close-up of the 12 C mass fraction after 10 
ms of evolution of the cold model at 0.5 (a), 2 (b), 4 (c), and 7.5 cm zone -1 (d) resolutions. 
The base of the burning layer is located in the bottom-most green region (just below the 
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magenta) in panel a. All four simulations give a well-mixed carbon region above the burning 
layer; the extent of the convective zone increases with decreasing resolution with the 7.5 cm 
zone -1 simulation's convective zone extending 30% further than the 0.5 cm zone -1 simula- 
tion's convective zone. The amount of convective undershoot — the tendency of material to 
penetrate below the burning layer — is much more sensitive to resolution. The 0.5 cm zone -1 
simulation shows very little evidence of undershooting while the 7.5 cm zone -1 simulation 
has an undershoot region that is larger in physical extent than its corresponding convective 
region above the burning layer. For all of the studies described below, we use a resolution 
of 0.5 cm zone -1 in the burning layer. 



4.2. Effects of Thermal Diffusion on the Burning Layer 

As explained in Section [TJ the burning front during an XRB likely propagates as a 
subsonic flame, the speed of which is regulated by the rate of thermal diffusion across the 
front. At the resolution required to resolve the thin burning layer (see previous section) it is 
currently intractable to evolve the system until flame ignition. We can, however, investigate 
the effects of thermal diffusion on the stable burning that occurs in the burning layer. Here we 
focus on the hot model instead of the cold model because it has the larger thermal gradient— 
and hence diffusive heat flux — at the base of the accreted layer. For this simulation, we use 



the new adaptive mesh refinement (AMR) capability in MAESTRO (iNonaka et al.ll2010l ). using 
2 levels of refinement and ensuring that the entire convective region is at the finest level 
of refinement with resolution 0.5 cm zone -1 . Figure |5] shows these effects in (-Hnuc) max 
(solid lines) and the location of this maximum (dashed lines) as a function of time at early 
times both with (green) and without (blue) thermal diffusion. We note that the location of 
(-^nuc) max is always at the finest level of refinement. The {H TLUC ) max evolution is similar for 
both cases with the magnitude in general being slightly larger for the case of no diffusion. 
The initial spike in {H nnc ) max at t ~ 0.25 ms is due to the fact that, initially, there is 
no established fluid flow that can advect away the energy released from nuclear reactions 
(see the discussion in Section I4T5]) . Over the next 3 ms, the location of (H nuc ) max for both 
simulations moves radially inward at a rate of ~ 2.9 x 10 3 cm s . Around t = 3.25 ms, the 
inward radial progression of the location of (-^nuc) max for the simulation with no diffusion 
significantly slows to ~ 900 cm s -1 . For the remainder of the simulation, the case with 
thermal diffusion shows no such slowdown — heat transported radially inward via diffusion 
expands the lower boundary of the convective zone, which mixes fresh fuel to slightly deeper 
layers. By the end of the simulations, the case that included diffusion had an (-Hnuc) max 
that occurred ~ 4 cm deeper within the atmosphere than in the case without diffusion. It 
should be noted that the typical standard deviation in the location of (-£^nuc) max for the case 
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without diffusion is of order 2 cm; this suggests that perhaps thermal diffusion plays a role 
in regulating the location of maximum nuclear burning, but further evolution is needed to 
make statistically significant claims. 



4.3. Expansion of Base State due to Heating 

Having a dynamical base state allows us to capture the large-scale expansi on of the 



atmos phere due to heating from nuclear reactions. This differs from the work by iLin et al. 



( 120061 ). which had a time- independent base state and did not model the top of the accreted 
atmosphere due to numerical complications with their algorithm. Figure E] shows the ratio of 
the base state density to that of the initial (t = 0) base state density profile near the surface 
of the atmosphere for the cold model. We define the surface to be where po = Pcutofr- The 
vertical dashed lines represent the location of the surface for each time- value. After 26.6 ms of 
evolution, the base state has responded to heating from nuclear reactions approximately 4.5 
m below the surface by expanding 3.5 cm. The lower Mach number flow in the cold model 
compared to the hot model allows for longer-term evolution of the system and therefore 
larger expansion of the atmosphere. 

The extent of the expansion is rather small at these early times. However, as the system 
progresses towards outburst the energy generation and, therefore, the rate of expansion in- 
creases. As the system expands, the po profile changes, which can affect the dynamics in the 
convective region. Additionally, as the atmosphere expands, the burning layer becomes less 
degenerate, which may be important for the nucleosynthesis during the outburst. Further- 
more, a proper modeling of this expansion during the peak of a PRE burst model may help 
pinpoint the location of the photosphere with respect to the stellar radius at touchdown, a 
quantity that plays an i mportant role in usin g XRBs to measure the mass and radius of the 



underlying neutron star ( jSteiner et al.ll2010l . for example) 



4.4. Effects of the Volume Discrepancy Term 

In Section 13.11 we explained that the thermodynamic pressure may drift from the base 
state pressure. To correct for this drift, we introduced the volume discrepancy term in 
equation ffTB"]) , which drives the thermodynamic pressure towards the base state pressure. 
We focus our attention here on the hot model because it shows a more dramatic drift 
of the thermodynamic pressure from the base state pressure. Figure [7] shows the volume 
discrepancy term in action by examining the percent difference between the base state and 
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thermodynamic pressures as a function of time for various values of / for the hot model. 
The top panel shows the maximum value whereas the bottom panel shows the average value 
of this percent difference; both the peak and average values show the same trend for a given 
value of /. After the initial adjustment of the system, the average drift for the case of no 
volume discrepancy correction (/ = 0) increases approximately linearly at ~ 0.1% per ms of 
evolution. Including the correction term restricts the temporal- and spatial-averaged value 
of the drift to < 0.02%. 

For nonzero /, the oscillatory behavior in the drift is due to the fact that the system 
may slightly over-correct the thermodynamic pressure in a given time step and then recover 
in the next step. A larger value of / causes a stronger driving of the drift, which tends to be 
more oscillatory. In addition, a larger value of / appears to be correlated with larger spikes 
in the drift. The top panel of Figure [8] shows a closeup of the 0(1) error seen in the / = 0.3 
curve in Figure [7J The location of the maximum drift is also plotted in the top panel; the 
large spike in the drift occurs just below the burning layer at r = 366.25 cm. The bottom 
panel of Figure [8] shows the corresponding maximum energy generation rate, which also 
contains a spike that is coincident with the spike in the drift — a large deposit of energy on a 
short timescale causes the thermodynamic pressure to get out of sync with the hydrostatic 
base state pressure. The increase in the energy generation rate is due to a fluid parcel rich 
in 4 He fuel being brought into a region of high temperature via the turbulent convection. 
The duration of this transient behaviour is very short: 9 time steps or ~ 6.4 x 10 -7 s. 
The selection of an appropriate non-zero value for / is a problem-specific endeavor, but the 
chosen value has little effect on the dynamics of the convective flow field. For the simulations 
presented below we use a volume discrepancy correction value of / = 0.3, which is based 
on the results of several test runs and past experience with comparing the results to the 
/ = case. We will continue to study if and how the chosen value of / affects the long term 
development of the convective field for this specific problem. 



4.5. Convective Dynamics 

The adiabatic excess, AV, — with 

AV = V - V s , (19) 

where the actual thermal gradient is 

d In T/dr 
d In p/dr 

and the adiabatic thermal gradient is V s = (d In T/d In p) s with the subscript s meaning along 
an adiabat — is used to gauge the evolution of the convective zone for the Schwarzschild 
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instability criterion. Under this criterion, a fluid element is unstable to thermally driven 
convection when AV > and is stable for AV < 0. The first term in (fT9l) is calculated 
using finite differences of the temperature and pressure profiles along the radial direction. 
The second term in f[T9"j) depends solely on the ther modynamics of the equation of state. It 



is related to the second adiabatic exponent, T 2 (see I Cox &: Giulil (119681 ) Chapter 9): 



dlnp 

All three of the adiabatic exponents are related: 

ri r 2 



(20) 



r 3 -i r 2 -i 



(21) 



where T 3 — 1 = (dhiT/dln p) s and Ti was defined in Section T3.ll Writing the equation of 
state as p = p(p, T) and expanding the differential dp, we find the relation 

T ,-l = ±^> (22) 
Tpr/P 

along an adiabat. Our equation of state only returns r 1; but combining this with (121 j) and 
allows us to solve for the adiabatic thermal gradient and hence the adiabatic excess. 



Figure [9] shows the early evolution of AV for the cold model. Each plot covers the 
spatial range (0 < x < 256 cm, 350 cm < r < 700 cm) to focus on the convective region. 



The stripes in the initial conditions, Figure 9(a), are due to small interpolation errors from 



mapping the initial data onto the two-dimensional grid. The initial adjustment of the system 



seen in Figure 9(b) causes a mixing of stable (blue) and unstable (red) fluid elements. This 
transient adjustment phase occurs for two reasons: 1) the initial conditions were based on 
a parameterization of convection in one dimension and the system now needs to adjust to a 
two-dimensional convective zone, and 2) the initial perturbation does not have an established 
convection zone and the system needs a short amount of time to build up a flow pattern 
associated with the perturbation. This results in mixing that produces a region that is 
marginally convective (AV ~ 0; white) with localized pockets of stable and unstable fluid 



elements as seen in Figure 9(c), At later times, these pockets further localize into vortices 
whose circulation gives rise to roughly circular regions of nonzero adiabatic excess — with 
one hemisphere that is stable and the other which is unstable — that are advected with the 
flow before dispersing into the ambient medium on subconvective timescales, ~ 10~ 4 s. The 
vortices are always associated with an adiabatic excess pattern that has an unstable (red) 
bottom and a stable (blue) top unless two vortices are merging and interacting, in which 
case the stability distribution becomes skewed. 



- 18 - 



Figure [TO] shows AV for the same simulation as in Figure M but at later times. The 
boxes in these plots outline a single long-lived vortex that forms around t = 18.5 ms, Figure 
10(a), and lasts throughout the remainder of the simulation. Formation of this vortex is 



correlated with the formation of stronger filamentary structures, which are most clear in 



Figures 10(d), 10(e) and 10(f). These filaments appear to wrap around the solitary vortex 
and restrict the main formation of smaller vortices to the lower boundary of the convective 
region. 

Another way to quantify the convective region is to look at the ratio V/V s . From 
(fT9]) we see that the system is unstable to convection under the Schwarzschild criterion 
when V/V s > 1. The Schwarzschild criterion, however, does not consider the effects of 
composition gradients that may help stabilize the material against convection; for this we 
need to consider the Ledoux criterion for instability 



V - V L > 0, 



(23) 



where the Ledoux thermal gradient is (see, for example, iKippenhahn fc Weigertl (119941 )) 

dlnXi/dr 



v.-E 



d In p/dr 



and the second term above is evaluated via finite differences of the composition and pressure 
profiles. As with the Schwarzschild criterion, one can look at the ratio V/Vl, which is 



greater than unity if the material is unstable to Ledoux convection. Figure 11(a) shows the 
above ratios for the average thermal gradients for the initial configuration (left) and after 
the system has evolved for t = 23.5 ms (right); the black line is for the case of Schwarzschild 
criterion convection, while the red line is for Ledoux convection. The dashed horizontal line 
marks the boundary for stability against convection. Where the curves lie above this line, 
the configuration is unstable; when convection is efficient, the curves should lie very near 
the horizontal line. For both the initial condition and at late times, the Schwarzschild curve 
and the Ledoux curve are well matched except near the edges of the convective region where 
composition gradients cause the two curves to deviate slightly. This is most noticeable in 
the initial configuration at the upper boundary where there is a sharp jump in composition, 
which was not smoothed (see Figured]). Of interest in the plot at t = 23.5 ms is the feature 
at r = 450 cm, which has an unstable bottom and a stable top; this is consistent with the 
vortices in Figures M and [TU] which had red bottoms and blue tops. We define the edge 
of the convective region to be where V/V s , V/Vl = 0.75. This particular value of 0.75 
was chosen to be sufficiently small enough to rule out false positives from strong pockets of 
stability from, for example, vortices within the convective region, but also large enough to 



rule out any fluctuations at the boundaries due to overshoot. Figure 11(b) shows in grey 
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the extent of the convective region as a function of time with respect to the full domain for 
both the Schwarzschild (left) and Ledoux (right) instability criteria. The horizontal dashed 
lines mark the initial location of the lower and upper boundaries. The overall expansion 
of the upper boundary for the Schwarzschild (Ledoux) criterion is 32.0 (29.5) cm in 30 ms 
of evolution; the lower boundary expands downward by 9.5 (6.0) cm in the same time. At 
late times, the upper boundary of the convective region has a much smoother composition 
transition than the lower boundary, therefore, the Schwarzschild and Ledoux criteria are 
much better matched at the upper boundary than the lower. Nevertheless, both the Ledoux 
and Schwarzschild criteria yield similar results when used to determine the extent of the 
convective region. In terms of column-depth, the convective zone after 30 ms of evolution 
spans the region 2.2 x 10 7 g cm -2 < y < 2.6 x 10 8 g cm -2 for both instability criteria. 

For comparison, Figures [T2] and [13] show the 12 C mass fraction with velocity vectors 
for the same simulation and at the same times as in Figures [9] and [TUJ, respectively. These 
figures clearly show the association of vortices with the circular regions of nonzero adiabatic 
excess seen in Figures [9] and [101 The initial adjustment of the system causes mixing that 
smooths the slight over- abundance of 12 C at the base of the accreted layer present in the 
initial model (see Figured]). At late times, the convective region is very well-mixed, and 
the 12 C mass fraction is nearly laterally homogeneous. Furthermore, the circulation pattern 
associated with the long-lived vortex outlined in Figure [10] has grown to a large fraction 
of the convective zone and is self-interacting because of the periodic boundary conditions. 
The tendency of the system to form a single dominant vortex from smaller vortices is a 
feature of two-dimensional simulations. In three dimensions, the turbulent energy cascade 
moves from large to small scales; large vortical structures break down into smaller structures 
that are eventually dissipated by viscous effects. In two dimensions, as is the case here, 
the turbulent energy cascade is reversed — small vortical structures merge together to form 
a single dominant vortex. 

The circulation is counter-clockwise for the large, long-lived vortex; this causes a region 
with positive x-velocity below and a region of negative x-velocity above the vortex center. 
The positive x-velocity region extends all the way to the lower convective boundary where 
it causes shearing of the 4 He/ 12 C-rich region with the underlying 56 Fe region. Consequently, 
Figure [TH shows that some of the underlying 56 Fe neutron star material is churned up into 
the convective region where it is mixed with the rest of the convective material. The left 
panel shows average 56 Fe mass fraction profiles starting with the initial model abundance 
(thick solid line) through t = 30 ms (thick dashed line); the intermediate thin solid lines 
show profiles at the times used in Figures M and [TU] By t = 5 ms, the 56 Fe is fairly well- mixed 
in the convective region. The right panel shows the total mass of 56 Fe in the region defined 
by the initial convective zone. The greatest growth in the total mass occurs, as expected, in 
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the initial adjustment (t < 0.6 ms) and then flattens until large enough structures form such 
that there is sufficient shearing occurring at the base of the convective boundary. There is 
only a slight increase in the growth rate for the 56 Fe mass around t = 18.5 ms where the 
long-lived vortex first appears. This is due to the fact, mentioned above, that as the system 
evolves it goes from many small vortices to a few large, dominant vortices. It is only when 
the circulation pattern of a particular vortex is large enough to strongly interact with the 
lower convection boundary that we get the shearing and enrichment of the convective region; 
this occurs around t ~ 5 ms. The addition of 56 Fe to the convective region has a small but 
noticeable effect on the conductivity; for example, a displacement of ~ 1% 4 He for 56 Fe near 
the base of the accreted layer, with all other things being constant, gives a ~ 4% decrease in 
conductivity. This could play an important role in adjusting the flame speed once ignited. 

Figure \T5\ shows the evolution of the maximum value of H nnc throughout the duration 
of the cold model simulation. The inset plot shows the early adjustment phase mentioned 
above. The initial jump in H nnc is due to the fact that there is no well established flow field 
that can efficiently advect away the energy released from reactions. Once the convective zone 
is well established, the energy generation rate relaxes before making its steady climb. We 
note that we have not yet achieved runaway — the rise in energy generation rate is still linear. 
This climb is temporarily interrupted by a couple of spikes similar to those seen in Figure 
[8] when fresh fuel is advected to a relatively hot region and burned quickly. Although well 
organized at later times, the convective fluid flow is slow with respect to the sound speed. 
Figure [TB] shows the maximum Mach number in the computational domain as a function of 
time; this value never exceeds 0.08 in our simulation. The average value of the Mach number 
in the convective region rarely exceeds ~ 0.02 during our 30 ms simulation. 



5. Conclusions 

We have described some of the challenges and important concepts to keep in mind when 
performing multidimensional simulations of XRBs. The major results can be summarized as 
follows: 

• To get a system that is much closer to thermal instability in multiple dimensions, the 
semi-analytic one-dimensional models should augment the local cooling rate estimate, 
01]), to include cooling due to convection. 

• Properly resolving the burning layer using the initial models considered here requires a 
spatial resolution of 0.5 cm zone -1 , which is an order of magni tude higher tha n what has 



been presented in the literature for multidimensional models ( iLin et al.ll2006l ). It should 
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be noted that our initial models differ in the underlying neutron star's composition — 
their 12 C opposed to our 56 Fe — and their models were evolved in one dimension through 
several bursts before being mapped into multiple dimensions. 

• Under-resolving the burning layer leads to dramatic convective undershoot and the 
burning tends to die out. 

• At the early times simulated here, the inclusion of thermal diffusion has little effect on 
the maximum energy generation rate, but does perhaps affect the depth at which this 
maximum occurs. 

• The MAESTRO algorithm we use allows us to capture the expansion of the atmosphere 
due to heating, which will be important in the modeling of PRE burst sources. 

• The average thermal gradient in the convective region is nearly adiabatic but there are 
localized pockets and filamentary structures that are either super- or sub-adiabatic. 

• The strong convection interacts with and churns up the underlying neutron star ma- 
terial, which slightly alters the conductivity of the convective region. 

The initial selection of a value to use for the volume discrepancy term in our simulations 
was based on experience with other applications. As we showed in Section 14.41 the value 
used for the long duration simulation in this paper, / = 0.3, may not the optimal choice 
for the XRB problem. Further investigation is required to determine which factors affect 
the appropriate value of /, and to determine if the spikes in the drift of the thermodynamic 
pressure from the base state pressure are simply numerical artifacts of a poor choice of /. 

The width of the computational domain used in our simulations is adequate for the 
early evolution of the system; the size of any individual convective cell is initially small with 
respect to the width of the domain. As the system evolves and the convection becomes more 
established, the cells grow in size. The nature of vorticity in two dimensions is such that the 
smaller vortices merge to form a single vortex. In our simulations the cells grow to become a 
significant fraction of the domain width and the flow becomes dominated by a single vortex 
that interacts with itself through the periodic boundary conditions. By selecting a wider 
computational grid, we could delay the formation of a single, dominant vortex. Ideally the 
computational domain should be several pressure scale-heights wide so that we should form 
multiples of these convective cells that dominate the flow for an extended period of time 
before merging into a single vortex. Given our strict resolution requirements, such a setup 
was computationally infeasible. 
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We plan to further investigate some of the topics mentioned above in future work while 
studying mixed H/He bursts. In such bursts the majority of the energy release comes from 
burning hydrogen; the nuclear reaction rates involved in such burning are less temperature 
sensitive than the 3- a rate used in the current paper. This may allow for a relaxed resolution 
requirement for properly resolving the burning layer because the energy generation rate 
profile should not be as sharply peaked as we have seen in our studies. This would allow 
for longer time evolution, which may allow us to say something about whether or not the 
convective zone extends all the way to the photosphere near outburst. We will also be able 
to simulate larger domains where we could address the effects of domain size on the long- 
term evolution of the convective region and its 56 Fe enrichment. Furthermore, we will begin 
investigating the effects of unprecedented three-dimensional simulations of the convection 
that precedes the outburst in an XRB and compare its properties to our two-dimensional 
studies. All of these simulations will rely heavily on MAESTRO's new AMR capability to 
reduce computational cost for long-term evolution. We will also investigate the effects of 
using different methods of c alcula ting opacities as well as updated reaction rates, such as 
the 3-a rate of lOgata et all (120091 ). 
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Fig. 1. — Kepler-supplemented cold (solid lines) and hot (dashed lines) models as described 
in the text. Energy release from nuclear burning at the base of the 4 He layer has caused the 
temperature to rise. The cold model is evolved to a peak Tb ase = 3.67 x 10 8 K and the hot 
model is evolved to a peak Tb aS e = 5.39 x 10 8 K. The black vertical lines indicate the location 
of the anelastic cutoff while the grey vertical lines indicate the location of the beginning of 
our sponge forcing term for each of the models (see Section I3TT1 ) 
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Fig. 2. — Initial sponge profile for the cold model where r sp = 680 cm and r tp = 844 cm. 
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Fig. 3. — Average of H nnc as a function of height for the cold model at various resolutions at 
t — 1 ms. Note that the vertical axes of the inset plots are in a logarithmic scale. For c larity , 
the top panel shows simulations which use the same resolutions as in the iLin et al.l (120061 ) 
study and the bottom panel shows more resolved simulations. The peak of the profile at 
0.5 cm zone -1 resolution is qualitatively similar to the peak of the profile at 0.25 cm zone -1 
resolution. 
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Fig. 4. — Effects of under-resolving convection for the cold model. Plotted is the 12 C mass 
fraction after 10 ms of evolution for various resolutions: a) 0.5, b) 2, c) 4, and d) 7.5 cm 
zone -1 . Each figure shows the same region of physical space and has dimensions 256 cm 
x 1024 cm. The coarse resolution simulations show an extended convective zone and a 
significant amount of convective undershoot. 
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Fig. 5. — Evolution of (H nuc ) max (solid lines) and its vertical location (dashed lines) as a 
function of time for the hot model both with (green) and without (blue) thermal diffusion. 
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Fig. 6. — Expansion of the base state due to heating. Plotted is the ratio of base state 
density, p , to the initial (t = 0) base state density, po.mit, near the surface of the atmosphere 
for the cold model. We define the surface to be where p = p cuto fi and it is represented by 
the vertical dashed lines. The base state has expanded 3.5 cm in 26.6 ms of evolution. 
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Fig. 7. — Effects of the volume discrepancy factor as characterized by the percent difference 
between the thermodynamic pressure as given by the equation of state, Peos ; an d the base 
state pressure, po, for the hot model. The top panel shows the maximum value whereas the 
bottom panel shows the average value of the percent difference in the computational domain. 
Note the different vertical scales between the two plots. 
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Fig. 8. — Closeup of the 0(1) spike in the maximum value of the / = 0.3 drift as seen 
in the top panel of Figure [7J The top panel shows the drift value and its location in the 
domain; the bottom panel shows the maximum energy generation rate. The large amount 
of energy released from the burning spike causes the thermodynamic pressure to differ from 
the hydrostatic base state pressure and therefore a spike in the drift. 
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(d)t = 5 ms (e) t = 7.5 ms (f) t = 10 ms 



Fig. 9. — Colormap plot of the evolution of the adiabatic excess, AV, in the convective 
region for the cold model. 
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(d) t = 25 ms (e) t = 26 ms (f) t = 28 ms 



Fig. 10. — Same as Figure M but at later times. The boxes show the location of a single 
feature that, once formed, lasts for the remainder of the simulation. 
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(b) Convective Region Extent 

Fig. 11. — Analysis of the extent of the convective region. Panel (a) shows the convective 
profiles for both the Schwarzschild and Ledoux instability criteria at two different times. 
Panel (b) shows the extent of the convective region as a function of time as determined by 
both instability criteria. 




Fig. 12. — Colormap plot of 12 C mass fraction with velocity vectors for the same region and 
times as shown in Figure [91 
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Fig. 13. — Colormap plot of 12 C mass fraction with velocity vectors for the same region and 
times as shown in Figure [10J 
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Fig. 14. — Plots showing the Fe enrichment of the convective region. The left panel shows 
the evolution of the average 56 Fe mass fraction starting from the initial model distribution 
(solid thick line) and ending after 30 ms of evolution (dashed line); the thin grey lines show 
the evolution at the intermediate times shown in Figures M and [10J The right panel shows 
the total mass of 56 Fe in the convective region as a function of time. Note the log scale of 
the horizontal axis in the right plot. 
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Fig. 15. — Plot of the maximum H nnc in the cold model simulation as a function of time. 



The inset plot shows the early adjustment phase associated with Figures 9(b) and 9(c) The 
spikes are similar to that seen in the bottom panel of Figure [SJ and are caused by the rapid 
burning of fresh fuel as it is brought into the burning layer by the turbulent convection. 
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Fig. 16. — Plot of the maximum Mach number in the cold model simulation as a function 
of time. The slow convective flow justifies the use of a low Mach number approximation 
method. 
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A. Thermal Diffusion 



Here we describe the changes from Paper V due to the inclusion of the thermal conduc- 
tion term in equation (Q. The boldface notation refers to specific steps in the algorithm, 
which are fully described in Appendix A. 4 in Paper V. 

Applying the chain rule to the equation of state, h = h(po,T,X), we note that the 
temperature gradient can be expressed as 



Whenever we require an explicit computation of the thermal conduction term, we use this 
formulation. In the edge state prediction (Steps 4H and 8H), we add an explicit con- 
tribution of the thermal conduction term to the forcing. Also, whenever we compute the 
expansion term S, we include the thermal conduction contribution. To compute thermody- 
namic derivatives of this term, we use h, X, and po as inputs to the equation of state. We 
account for thermal diffusion in the cell update step by replacing Steps 41 and 81 with the 
semi-implicit approach described below. 

Step 41. Diffuse the enthalpy through a time interval of At. 

Compute k[^,Cp \ and £jp from p^\T^\ and xjp as inputs to the equation of state. 
We denote the result for enthalpy in Step 4H of Paper V as (ph)^ 2 '^'* rather than (ph)^'* 
to indicate that we are about to account for thermal diffusion and define 

p (2'),* = p (2),*. The 

update is given by 




(Al) 




(A2) 
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which is numerically implemented as a diffusion equation for h^'*, 

At k (1) \ At h [l) 

c p J c p 

_^V.(«V P f.* + «vrf'), (A3) 

\ C p Cp J 

Then, update temperature using the equation of state: T^'* = T (p( 2) '*, /i (2) '*, X^'*) . 

Step 81. Diffuse the enthalpy through a time interval of At. 

Compute kfa'*, Cp 2 ' 1 '*, and £^ , from p^'*, T^'*, and X^'* as inputs to the equation of 
state. We also denote the result for enthalpy in Step 8H of Paper V as (ph)^ 2 '^ rather than 
(ph)@) to indicate that we are about to account for thermal diffusion and define p^ 2 '^ = p( 2 \ 
The update is given by 

At I V 2) '* k {1) 



.(2),*,(2) I * t(l) t (l) 



-f v ' t"1 2) + ^v p < "i 

\ Cp Cp 



which is numerically implemented as a diffusion equation for /i^ 2 \ 
A/ £* {2) '* \ A/ k {1) 

p (2) -f v -%^vU (2) = W (2,) + f v • %v^) 



A/ / .(2),*, (2),* /^M 1 ) 

fc \ C P °v 

,(2),* 1.(2),* /^M 1 ) 



^ Vp «, + ^ Vp a),, (A5) 



Then, update the temperature using the equation of state: T (2) = T (p (2) , h (2 \ xf ^ 
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A.l. Diffusion Solver Test 

This problem is designed to test the accuracy of our implementation of an implicit 
solver for the diffusion of a two-dimensional Gaussian enthalpy pulse. That is, we are only 
concerned with the diffusive term in (19]): 



d( P h) 
dt 



V • (fct h VT) • 



(A6) 



To easily compare with an analytic solution (see, for example ISwesty fc Myral (120091 ) for an 
analogous example for a radiation-hydrodynamics code) we assume the thermal conductivity 
to be constant: fcth = 10 7 erg K cm -1 s _1 . Note that this does not fully test the predictor- 
corrector aspect of the method outlined in Appendix |X] because in this simplified problem 



fcW'* = We also assume an ideal gas with X(He 4 ) = 0.5, X(C 12 ) = AT(Fe 5B ) = 0.25 and 
ratio of specific heats 7 = 5/3. Furthermore, we are not concerned with any hydrodynamic 



motions so we keep the density fixed. We can then express (1A6I) in a simpler form 



dh 
dt 



DV 2 K 



(AT) 



where D = k t h/ (p c P ) is the diffusion coefficient. 

Given the initial conditions for the two-dimensional pulse, 



h(r, t = to) — (^p ~~ ^0) x ex P 



r o 



where h p , ho, r = (xo,yo), and to are the peak enthalpy, ambient enthalpy, location of the 
center of the peak, and time from which the system has evolved respectively, the analytic 
solution takes on the form 



h(r,t) = (h p - ho) 



( <0 ) 


exp 


(- 


r - r 


2 \ 


\t + t ) 


1 AD (t + t ) I 



h , 



(A9) 



where t is the evolved time. 



We solve this problem on a Cartesian grid of size 4 cm x 4 cm with the following 
parameters: h p = 10.0 erg g -1 , h = 1.0 erg g" 1 , r = (2.0 cm, 2.0 cm), t = 0.1 s, and 
p = 1.0 g cm -3 . For the density and composition used in this test, we obtain a diffusion 
coefficient of D = 0.32 cm 2 s _1 . As motivated in Appendix |Aj our implicit solve uses a 
Crank-Nicholson scheme that is second order accurate in space and time. 

Figure [T7] shows an example of the initial enthalpy pulse and its evolution through 
t = 0.4 s on a 1024 x 1024 grid with fixed time step At = 10~ 3 s. Note that as the pulse 
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expands it begins to interact with the edges of the computational domain and the symmetry 
of the Gaussian peak is broken. Figure [18] shows the computed average enthalpy as a function 
of radius (X's) compared to the analytic solution (lines) for the same test problem shown in 
Figure [T71 Again, excepting boundary effects the numerical and analytic solutions are well 
matched. 

To check the convergence of the algorithm we ran simulations with various resolutions 
and compared the errors. To measure the error in the simulation, we use the L\ norm of 
the difference between the analytic and numeric solutions normalized to the L\ norm of the 
analytic solution, which we define as e: 

~ \\h(r,t-)\\ Ll /'(>;.,•/"') ' 1 J 

where h(rij,t m ) is the analytic solution at ry = ((a?j — xo) 2 + (yj — yo)) 1 ^ 2 and time t m = 
mAt and h™^ is the numeric solution at (xi,yj) and time t m . We further define the conver- 
gence rate, a, by comparing the value of e at the current resolution to the value of e at a 
finer resolution simulation: 

a = \og 2 [-^-\ (All) 

V I 6 \ finer/ 

For our comparisons, we take "finer" to mean a simulation with twice the resolution; to 
compare the simulations at the same physical time, the finer simulation must have evolved 
through twice the number of time steps as the coarser simulation. If our algorithm truly is 
second order accurate in space and time then a should equal 2. Table [Q shows the values 
of e and the convergence rate for various resolutions at t = 0.08 s; for a, the norm in the 
current column is compared to the norm of the finer resolution simulation in the column to 
its right. Our values of a agree very well with the expected value. 



Table 1: Reduced Li norms and convergence rate for the diffusion test problem at t = 0.08 
s. 



128 x 128 Error 


256 x 256 Error 


512 x 512 Error 


1024 x 1024 Error 


At = 0.008 s 


At = 0.004 s 


At = 0.002 s 


At = 0.001 s 


e 8.64 x 10" 5 


2.16 x 10~ 5 


5.39 x 10" 6 


1.35 x 10~ 6 


a 2.0012 


1.9999 


1.9988 
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Fig. 17. — Time evolution of the diffusion of a two-dimensional Gaussian pulse of enthalpy 
as described in the text. The value of time displayed is the evolution time, t. This simulation 
was run with a 1024 x 1024 grid with time step size At = 0.001 s. Excepting edge effects 
near the domain boundary, the numerical solution maintains its axisymmetric form about 
the center of the pulse at (x,y) = (2.0,2.0). 
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Fig. 18. — The average of enthalpy as a function of radius from the center, (x, y) = (2.0, 2.0), 
of a two-dimensional Gaussian pulse. The X's are data from the numerical solution at the 
shown times. The lines represent the analytic solutions as given by f ]A9j) . The numerical 
solution tracks the analytic solution very well except when the pulse has diffused enough 
that it begins to interact with the boundaries of the computational domain as seen in the 
inset plot. 



